home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpkdrv.zip
/
SS.FOR
< prev
next >
Wrap
Text File
|
1984-01-07
|
19KB
|
656 lines
C MAIN PROGRAM
INTEGER LUNIT
C ALLOW 5000 UNDERFLOWS.
C CALL TRAPS(0,0,5001,0,0)
C
C OUTPUT UNIT NUMBER
C
LUNIT = 6
C
CALL SSITS(LUNIT)
STOP
END
SUBROUTINE SSITS(LUNIT)
C LUNIT IS THE OUTPUT UNIT NUMBER
C
C TESTS
C SSICO,SSIFA,SSISL,SSIDI,SSPCO,SSPFA,SSPSL,SSPDI
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK SSICO,SSISL,SSIDI,SSPCO,SSPSL,SSPDI
C EXTERNAL SSIXX,SMACH
C BLAS SAXPY,SSWAP,SASUM,SCOPY
C FORTRAN ABS,AMAX1,FLOAT,IABS
C
C INTERNAL VARIABLES
C
REAL A(15,15),AINV(15,15),ASAVE(15,15),AJK,AJKP1
REAL AP(120),B(15),X(15),XEXACT(15)
REAL XP(15),T,Z(15)
REAL ANORM,AINORM,COND,COND1,DET(2),DETP(2)
REAL EN,ENORM,EPS,FNORM,ONEPX,Q(6),QS(6),RCOND
REAL RCONDP,RNORM,S,SASUM,SMACH,XNORM
INTEGER I,INERT(3),INERTP(3),IQ(6),J,JB
INTEGER K,KASE,KSING,KM1,KOUNT,KPFAIL,KPVT(15),KPVTP(15)
INTEGER KS,KSTEP,KSUSP(6),LDA,LUNIT,N,NPRINT
LOGICAL KPF
C
LDA = 15
C
C WRITE MATRIX AND SOLUTIONS IF N .LE. NPRINT
C
NPRINT = 3
C
WRITE (LUNIT,540)
WRITE (LUNIT,950)
C
DO 10 I = 1, 6
KSUSP(I) = 0
10 CONTINUE
KSING = 0
KPFAIL = 0
C
C SET EPS TO ROUNDING UNIT FOR REAL ARITHMETIC
C
EPS = SMACH(1)
WRITE (LUNIT,550) EPS
WRITE (LUNIT,530)
C
C START MAIN LOOP
C
KASE = 1
20 CONTINUE
C
C GENERATE TEST MATRIX
C
CALL SSIXX(A,LDA,N,KASE,LUNIT)
C
C N = 0 SIGNALS NO MORE TEST MATRICES
C
C ...EXIT
IF (N .LE. 0) GO TO 520
ANORM = 0.0E0
DO 30 J = 1, N
ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
30 CONTINUE
WRITE (LUNIT,700) ANORM
C
IF (N .GT. NPRINT) GO TO 50
WRITE (LUNIT,530)
DO 40 I = 1, N
WRITE (LUNIT,740) (A(I,J), J = 1, N)
40 CONTINUE
WRITE (LUNIT,530)
50 CONTINUE
C
C GENERATE EXACT SOLUTION
C
XEXACT(1) = 1.0E0
IF (N .GE. 2) XEXACT(2) = 0.0E0
IF (N .LE. 2) GO TO 70
DO 60 I = 3, N
XEXACT(I) = -XEXACT(I-2)
60 CONTINUE
70 CONTINUE
C
C SAVE MATRIX AND GENERATE R.H.S.
C
DO 90 I = 1, N
B(I) = 0.0E0
DO 80 J = 1, N
ASAVE(I,J) = A(I,J)
B(I) = B(I) + A(I,J)*XEXACT(J)
80 CONTINUE
X(I) = B(I)
XP(I) = X(I)
90 CONTINUE
C
C FACTOR AND ESTIMATE CONDITION
C
RCOND = -1.0E0
CALL SSICO(A,LDA,N,KPVT,RCOND,Z)
WRITE (LUNIT,820) (KPVT(I), I = 1, N)
C
C OUTPUT NULL VECTOR
C
IF (N .GT. NPRINT) GO TO 110
WRITE (LUNIT,760)
DO 100 I = 1, N
WRITE (LUNIT,770) Z(I)
100 CONTINUE
WRITE (LUNIT,530)
110 CONTINUE
C
C FACTOR PACKED FORM AND COMPARE
C
KPF = .FALSE.
K = 0
DO 130 J = 1, N
DO 120 I = 1, J
K = K + 1
AP(K) = ASAVE(I,J)
120 CONTINUE
130 CONTINUE
RCONDP = -1.0E0
CALL SSPCO(AP,N,KPVTP,RCONDP,Z)
IF (RCONDP .EQ. RCOND) GO TO 140
WRITE (LUNIT,840)
WRITE (LUNIT,880) RCOND,RCONDP
KPF = .TRUE.
140 CONTINUE
K = 0
KOUNT = 0
DO 160 J = 1, N
DO 150 I = 1, J
K = K + 1
IF (AP(K) .NE. A(I,J)) KOUNT = KOUNT + 1
150 CONTINUE
160 CONTINUE
IF (KOUNT .EQ. 0) GO TO 170
WRITE (LUNIT,840)
WRITE (LUNIT,890) KOUNT
KPF = .TRUE.
170 CONTINUE
C
C TEST FOR SINGULARITY
C
IF (RCOND .GT. 0.0E0) GO TO 180
WRITE (LUNIT,750) RCOND
KSING = KSING + 1
CALL SSIDI(A,LDA,N,KPVT,DET,INERT,Z,110)
WRITE (LUNIT,790) DET(1)
WRITE (LUNIT,810) INERT
GO TO 510
180 CONTINUE
COND = 1.0E0/RCOND
WRITE (LUNIT,570) COND
ONEPX = 1.0E0 + RCOND
IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,560)
C
C COMPUTE INVERSE, DETERMINANT, INERTIA, TRUE CONDITION
C
DO 200 J = 1, N
DO 190 I = 1, J
AINV(I,J) = A(I,J)
190 CONTINUE
200 CONTINUE
CALL SSIDI(AINV,LDA,N,KPVT,DET,INERT,Z,111)
AINORM = 0.0E0
DO 220 J = 1, N
DO 210 I = J, N
AINV(I,J) = AINV(J,I)
210 CONTINUE
AINORM = AMAX1(AINORM,SASUM(N,AINV(1,J),1))
220 CONTINUE
COND1 = ANORM*AINORM
WRITE (LUNIT,580) COND1
WRITE (LUNIT,790) DET(1)
WRITE (LUNIT,800) DET(2)
WRITE (LUNIT,810) INERT
C
C SOLVE A*X = B
C
CALL SSISL(A,LDA,N,KPVT,X)
C
IF (N .GT. NPRINT) GO TO 240
WRITE (LUNIT,590)
DO 230 I = 1, N
WRITE (LUNIT,780) X(I)
230 CONTINUE
WRITE (LUNIT,530)
240 CONTINUE
C
C MORE PACKED COMPARE
C
CALL SSPSL(AP,N,KPVTP,XP)
KOUNT = 0
DO 250 I = 1, N
IF (XP(I) .NE. X(I)) KOUNT = KOUNT + 1
250 CONTINUE
IF (KOUNT .EQ. 0) GO TO 260
WRITE (LUNIT,840)
WRITE (LUNIT,900) KOUNT
KPF = .TRUE.
260 CONTINUE
CALL SSPDI(AP,N,KPVTP,DETP,INERTP,Z,111)
IF (DETP(1) .EQ. DET(1) .AND. DETP(2) .EQ. DET(2))
* GO TO 270
WRITE (LUNIT,840)
WRITE (LUNIT,910) DETP
KPF = .TRUE.
270 CONTINUE
KOUNT = 0
K = 0
DO 290 J = 1, N
DO 280 I = 1, J
K = K + 1
IF (AP(K) .NE. AINV(I,J)) KOUNT = KOUNT + 1
280 CONTINUE
290 CONTINUE
IF (KOUNT .EQ. 0) GO TO 300
WRITE (LUNIT,840)
WRITE (LUNIT,920) KOUNT
KPF = .TRUE.
300 CONTINUE
C
C RECONSTRUCT A FROM FACTORS.
C
K = 1
310 IF (K .GT. N) GO TO 400
KM1 = K - 1
IF (KPVT(K) .LT. 0) GO TO 340
C
C 1 BY 1
C
IF (KM1 .LT. 1) GO TO 330
DO 320 JB = 1, KM1
J = K - JB
AJK = A(K,K)*A(J,K)
T = AJK
CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
A(J,K) = -AJK
320 CONTINUE
330 CONTINUE
KSTEP = 1
GO TO 370
340 CONTINUE
C
C 2 BY 2
C
IF (KM1 .LT. 1) GO TO 360
DO 350 JB = 1, KM1
J = K - JB
AJK = A(K,K)*A(J,K) + A(K,K+1)*A(J,K+1)
AJKP1 = A(K,K+1)*A(J,K) + A(K+1,K+1)*A(J,K+1)
T = AJK
CALL SAXPY(J,T,A(1,K),1,A(1,J),1)
T = AJKP1
CALL SAXPY(J,T,A(1,K+1),1,A(1,J),1)
A(J,K) = -AJK
A(J,K+1) = -AJKP1
350 CONTINUE
360 CONTINUE
KSTEP = 2
370 CONTINUE
C
C SWAP
C
KS = IABS(KPVT(K))
IF (KS .EQ. K) GO TO 390
CALL SSWAP(KS,A(1,KS),1,A(1,K),1)
DO 380 JB = KS, K
J = K + KS - JB
T = A(J,K)
A(J,K) = A(KS,J)
A(KS,J) = T
380 CONTINUE
IF (KSTEP .EQ. 2)
* CALL SSWAP(1,A(KS,K+1),1,A(K,K+